04_plot_obs_v_sim

library(dplyr)
library(data.table)
library(ggplot2)
library(GGally)
library(tidyr)
"%&%" = function(a,b) paste(a,b,sep="")
pre.dir <- "/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/"
my.dir <- pre.dir %&% "paper-reviewer-requests/OTD_simulations/"
obs.dir <- my.dir %&% "obs_otd_exp/"
sim.dir <- my.dir %&% "lmer.fits/"
obs.h2.dir <- pre.dir %&% "gtex-h2-estimates/"

See https://github.com/hwheeler01/GenArch/blob/master/paper-reviewer-requests/OTD_simulations/01_make_sim.r for how simulated expression phenotypes were generated.

compare h2 sim vs obs

tislist <- scan(my.dir %&% "ten.spaces.tissue.list", "c", sep='\n')
simlist <- scan(my.dir %&% "simlist", "c")

for(tis in tislist){
  h2obs <- read.table(obs.h2.dir %&% 'GTEx.' %&% tis %&% '.h2.all.models_FHSfdr0.05.Chr1-22_globalOtherChr.2015-10-06.txt', header=T,sep='\t') %>% select(tissue,ensid,gene,local.h2) %>% rename(obs_h2 = local.h2)
  for(i in 1:length(simlist)){
    sim <- simlist[i]
    h2sim <- read.table(my.dir %&% 'GTEx.' %&% tis %&% '_' %&% sim %&% '.h2.all.models_FHSfdr0.05.Chr1-22_globalOtherChr.2016-06-22.txt',header=T,sep='\t') %>% select(ensid,local.h2)
    colnames(h2sim) <- c('ensid','sim' %&% i)
    if(exists("h2join") == FALSE){
      h2join <- inner_join(h2obs,h2sim,by='ensid') 
    }else{
      h2join <- inner_join(h2join,h2sim,by='ensid')
    }
  }
  if(exists("allh2") == FALSE){
    allh2 <- h2join
  }else{
    allh2 <- rbind(allh2,h2join)
  }
  rm("h2join")
}


ggplot(allh2,aes(x=sim1,y=obs_h2)) + geom_point(alpha=1/4) + facet_wrap(~tissue,ncol=3) + theme_bw() + geom_smooth() + ggtitle(simlist[1]) + xlab("simulation h2") + ylab("observed h2")

ggplot(allh2,aes(x=sim2,y=obs_h2)) + geom_point(alpha=1/4) + facet_wrap(~tissue,ncol=3) + theme_bw() + geom_smooth() + ggtitle(simlist[2]) + xlab("simulation h2") + ylab("observed h2")

ggplot(allh2,aes(x=sim3,y=obs_h2)) + geom_point(alpha=1/4) + facet_wrap(~tissue,ncol=3) + theme_bw() + geom_smooth() + ggtitle(simlist[3]) + xlab("simulation h2") + ylab("observed h2")

ggplot(allh2,aes(x=sim4,y=obs_h2)) + geom_point(alpha=1/4) + facet_wrap(~tissue,ncol=3) + theme_bw() + geom_smooth() + ggtitle(simlist[4]) + xlab("simulation h2") + ylab("observed h2")

gallh2 <- gather(allh2,sim,sim_h2,5:8)
ggplot(gallh2,aes(x=sim_h2,y=obs_h2,color=sim)) + geom_smooth() + facet_wrap(~tissue,ncol=3) + theme_bw() + xlab("simulation h2") + ylab("observed h2")

ggplot(gallh2,aes(x=sim_h2,y=obs_h2,color=sim)) +geom_point(alpha=1/5) + geom_smooth() + facet_wrap(~tissue,ncol=3) + theme_bw() + xlab("simulation h2") + ylab("observed h2")

print(simlist)

## [1] "sim_exp_phenotype_errvar-ct_mult-1_seed-123"
## [2] "sim_exp_phenotype_errvar-ct_mult-2_seed-123"
## [3] "sim_exp_phenotype_errvar-ts_mult-1_seed-123"
## [4] "sim_exp_phenotype_errvar-ts_mult-2_seed-123"